
<!DOCTYPE html
  PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
   <head>
      <meta http-equiv="Content-Type" content="text/html; charset=utf-8">
   
      <!--
This HTML is auto-generated from an M-file.
To make changes, update the M-file and republish this document.
      -->
      <title>solar_array_3</title>
      <meta name="generator" content="MATLAB 7.8">
      <meta name="date" content="2011-12-10">
      <meta name="m-file" content="solar_array_3"><style type="text/css">

body {
  background-color: white;
  margin:10px;
}

h1 {
  color: #990000; 
  font-size: x-large;
}

h2 {
  color: #990000;
  font-size: medium;
}

/* Make the text shrink to fit narrow windows, but not stretch too far in 
wide windows. */ 
p,h1,h2,div.content div {
  max-width: 600px;
  /* Hack for IE6 */
  width: auto !important; width: 600px;
}

pre.codeinput {
  background: #EEEEEE;
  padding: 10px;
}
@media print {
  pre.codeinput {word-wrap:break-word; width:100%;}
} 

span.keyword {color: #0000FF}
span.comment {color: #228B22}
span.string {color: #A020F0}
span.untermstring {color: #B20000}
span.syscmd {color: #B28C00}

pre.codeoutput {
  color: #666666;
  padding: 10px;
}

pre.error {
  color: red;
}

p.footer {
  text-align: right;
  font-size: xx-small;
  font-weight: lighter;
  font-style: italic;
  color: gray;
}

  </style></head>
   <body>
      <div class="content"><pre class="codeinput"><span class="comment">%</span>
<span class="comment">% &#1057;&#1090;&#1088;&#1091;&#1082;&#1090;&#1091;&#1088;&#1072;, &#1089;&#1086;&#1076;&#1077;&#1088;&#1078;&#1072;&#1097;&#1072;&#1103; &#1087;&#1072;&#1088;&#1072;&#1084;&#1077;&#1090;&#1088;&#1099; &#1084;&#1077;&#1093;&#1072;&#1085;&#1080;&#1095;&#1077;&#1089;&#1082;&#1086;&#1081; &#1089;&#1080;&#1089;&#1090;&#1077;&#1084;&#1099;</span>
<span class="comment">%</span>
Model=struct;
<span class="comment">% &#1050;&#1086;&#1083;&#1080;&#1095;&#1077;&#1089;&#1090;&#1074;&#1086; &#1090;&#1077;&#1083;</span>
n=3;
Model.n=n;
<span class="comment">% &#1064;&#1072;&#1088;&#1085;&#1080;&#1088;&#1085;&#1099;&#1077; &#1074;&#1077;&#1082;&#1090;&#1086;&#1088;&#1099;</span>
<span class="comment">% &#1064;&#1072;&#1088;&#1085;&#1080;&#1088;&#1085;&#1099;&#1081; &#1074;&#1077;&#1082;&#1090;&#1086;&#1088; $c_{i,j}$ &#1085;&#1072;&#1095;&#1080;&#1085;&#1072;&#1077;&#1090;&#1089;&#1103; &#1074; &#1094;&#1077;&#1085;&#1090;&#1088;&#1077; &#1084;&#1072;&#1089;&#1089; &#1090;&#1077;&#1083;&#1072; &#1080; &#1079;&#1072;&#1082;&#1072;&#1085;&#1095;&#1080;&#1074;&#1072;&#1077;&#1090;&#1089;&#1103;</span>
<span class="comment">% &#1074; &#1096;&#1072;&#1088;&#1085;&#1080;&#1088;&#1077; $j$.</span>
Model.c=cell(n,n);
c=cell(n,n);
<span class="keyword">for</span> i=1:n
    <span class="keyword">for</span> j=1:n
        c{i,j}=[0;0];
    <span class="keyword">end</span>
    c{i,i}=[-0.5-0.005;(-1)^(i+1)*0.02];
    <span class="keyword">if</span> i&lt;n
        c{i,i+1}=[+0.5+0.005;(-1)^(i)*0.02];
    <span class="keyword">end</span>
<span class="keyword">end</span>
Model.c=c;
<span class="comment">% &#1052;&#1072;&#1090;&#1088;&#1080;&#1094;&#1099; &#1084;&#1072;&#1089;&#1089; &#1090;&#1077;&#1083; &#1089;&#1080;&#1089;&#1090;&#1077;&#1084;&#1099;. &#1056;&#1072;&#1089;&#1089;&#1084;&#1072;&#1090;&#1088;&#1080;&#1074;&#1072;&#1077;&#1084;&#1072;&#1103; &#1089;&#1080;&#1089;&#1090;&#1077;&#1084;&#1072; &#1087;&#1083;&#1086;&#1089;&#1082;&#1072;&#1103;, &#1087;&#1086;&#1101;&#1090;&#1086;&#1084;&#1091;</span>
<span class="comment">% &#1084;&#1072;&#1090;&#1088;&#1080;&#1094;&#1072; &#1084;&#1072;&#1089;&#1089; &#1082;&#1072;&#1078;&#1076;&#1086;&#1075;&#1086; &#1090;&#1077;&#1083;&#1072; &#1080;&#1084;&#1077;&#1077;&#1090; &#1074;&#1080;&#1076;:</span>
<span class="comment">% &lt;latex&gt;</span>
<span class="comment">% $$ M_i = \begin{pmatrix} m_i &amp; 0 &amp; 0 \\ 0 &amp; m_i &amp; 0 \\ 0 &amp; 0 &amp; J_{iz}</span>
<span class="comment">% \end{pmatrix}$$</span>
<span class="comment">% &lt;/latex&gt;</span>
<span class="comment">%</span>
Model.mass=cell(n,1);
mass=cell(n,1);
<span class="keyword">for</span> i=1:n
    mass{i}=[10 0 0; 0 10 0; 0 0 0.5];
<span class="keyword">end</span>
Model.mass=mass;

q0=zeros(n*2+3,1);
<span class="keyword">for</span> i=1:n
    q0(i)=(-1)^(i+1)*pi;
<span class="keyword">end</span>
q0(1)=pi/2;

[t,q]=ode113(@dqdtsol,[0 8],q0);
</pre><pre class="codeinput">dlmwrite(<span class="string">'result.csv'</span>, [t q(:,1:6)], <span class="string">','</span>);
</pre><pre class="codeinput">figure
hold <span class="string">on</span>
plot(t,q(:,1),<span class="string">'k-'</span>,<span class="string">'LineWidth'</span>,2);
plot(t,q(:,2),<span class="string">'k--'</span>,<span class="string">'LineWidth'</span>,2);
plot(t,q(:,3),<span class="string">'k-.'</span>,<span class="string">'LineWidth'</span>,2);
hold <span class="string">off</span>
legend(<span class="string">'\phi_1'</span>,<span class="string">'\phi_2'</span>,<span class="string">'\phi_3'</span>);
</pre><pre class="codeinput">figure
hold <span class="string">on</span>
plot(t,q(:,4),<span class="string">'k-'</span>,<span class="string">'LineWidth'</span>,2);
plot(t,q(:,5),<span class="string">'k--'</span>,<span class="string">'LineWidth'</span>,2);
plot(t,q(:,6),<span class="string">'k-.'</span>,<span class="string">'LineWidth'</span>,2);
hold <span class="string">off</span>
text(<span class="string">'interpreter'</span>,<span class="string">'latex'</span>);
legend(<span class="string">'\dot \phi_2'</span>,<span class="string">'\dot \phi_2'</span>,<span class="string">'\dot \phi_3'</span>);
</pre><pre class="codeinput">points=[0.005 1.005 1.005 0.005 0.005; -0.005 -0.005 -0.035 -0.035 -0.005];

M = moviein(size(t,1));
axis([-0.1 3.2 -1.2 2]);

hold <span class="string">on</span>;

p=cell(n,1);

maxY=0;

tmax=max(t);

<span class="keyword">for</span> i=1:tmax/500:tmax
    cla;
    a=0;
    dr=[0;0];
    <span class="keyword">for</span> j=1:3
        f=interp1(t,q(:,j),i);
        a=a+f;

        pp=getA(a)*[points(1,:); points(2,:)*(-1)^(j+1)];

        <span class="keyword">if</span> j&gt;1
            dr=dr+getA(a-f)*[1;-0.04*(-1)^j];
        <span class="keyword">end</span>

        line(pp(1,:) + dr(1), pp(2,:) + dr(2));
    <span class="keyword">end</span>

    f=getframe;

    <span class="comment">% M(:,i) = getframe;</span>
    <span class="comment">% imwrite(M(:,i),'img.png','PNG','Resolution',100);</span>
    <span class="comment">% mo= getframe;</span>
    <span class="comment">% mov=addframe(mov,mo);</span>
<span class="keyword">end</span>
</pre><p>Energy</p><pre class="codeinput">qabs=q*0;
v=zeros(size(q,1),2);
h=zeros(size(q,1),2);

T=0;
P=0;

<span class="keyword">for</span> i=1:n
    v=zeros(size(q,1),2);
    h=zeros(size(q,1),2);
    <span class="keyword">for</span> j=1:i
       qabs(:,i)=qabs(:,i)+q(:,j);
       qabs(:,n+i)=qabs(:,n+i)+q(:,n+j);
    <span class="keyword">end</span>
    <span class="keyword">for</span> j=1:i
       <span class="keyword">if</span> j==i
           b=-c{j,j};
       <span class="keyword">else</span>
           b=-c{j,j}+c{j,j+1};
       <span class="keyword">end</span>
       h=h+reshape(getA(qabs(:,j))*b,size(q,1),2);
       v=v-reshape(getA(qabs(:,j))*[0 1;-1 0]*b,size(q,1),2).*[qabs(:,n+j) qabs(:,n+j)];
    <span class="keyword">end</span>
    T=T+mass{i}(1,1)*(v(:,1).^2+v(:,2).^2)*0.5+mass{i}(3,3)*qabs(:,n+i).^2*0.5;
    P=P+Model.mass{i}(1,1)*9.81*h(:,2);
<span class="keyword">end</span>

figure;
E0=T+P;
plot(t,T+P-E0(1));
title(<span class="string">'Energy error'</span>);
</pre><p class="footer"><br>
            Published with MATLAB&reg; 7.8<br></p>
      </div>
      <!--
##### SOURCE BEGIN #####
%
% Структура, содержащая параметры механической системы
%
Model=struct;
% Количество тел
n=3;
Model.n=n;
% Шарнирные векторы
% Шарнирный вектор $c_{i,j}$ начинается в центре масс тела и заканчивается
% в шарнире $j$.
Model.c=cell(n,n);
c=cell(n,n);
for i=1:n
    for j=1:n
        c{i,j}=[0;0];        
    end
    c{i,i}=[-0.5-0.005;(-1)^(i+1)*0.02];
    if i<n
        c{i,i+1}=[+0.5+0.005;(-1)^(i)*0.02];
    end
end
Model.c=c;
% Матрицы масс тел системы. Рассматриваемая система плоская, поэтому
% матрица масс каждого тела имеет вид:
% <latex>
% $$ M_i = \begin{pmatrix} m_i & 0 & 0 \\ 0 & m_i & 0 \\ 0 & 0 & J_{iz}
% \end{pmatrix}$$
% </latex> 
% 
Model.mass=cell(n,1);
mass=cell(n,1);
for i=1:n
    mass{i}=[10 0 0; 0 10 0; 0 0 0.5];
end
Model.mass=mass;

q0=zeros(n*2+3,1);
for i=1:n
    q0(i)=(-1)^(i+1)*pi;
end
q0(1)=pi/2;

[t,q]=ode113(@dqdtsol,[0 8],q0);


%%

dlmwrite('result.csv', [t q(:,1:6)], ',');

%%
figure
hold on
plot(t,q(:,1),'k-','LineWidth',2);
plot(t,q(:,2),'kREPLACE_WITH_DASH_DASH','LineWidth',2);
plot(t,q(:,3),'k-.','LineWidth',2);
hold off
legend('\phi_1','\phi_2','\phi_3');

%%
figure
hold on
plot(t,q(:,4),'k-','LineWidth',2);
plot(t,q(:,5),'kREPLACE_WITH_DASH_DASH','LineWidth',2);
plot(t,q(:,6),'k-.','LineWidth',2);
hold off
text('interpreter','latex');
legend('\dot \phi_2','\dot \phi_2','\dot \phi_3');
%%

points=[0.005 1.005 1.005 0.005 0.005; -0.005 -0.005 -0.035 -0.035 -0.005];

M = moviein(size(t,1));
axis([-0.1 3.2 -1.2 2]);

hold on;

p=cell(n,1);

maxY=0;

tmax=max(t);

for i=1:tmax/500:tmax
    cla; 
    a=0;
    dr=[0;0];
    for j=1:3 
        f=interp1(t,q(:,j),i);
        a=a+f;
        
        pp=getA(a)*[points(1,:); points(2,:)*(-1)^(j+1)];
        
        if j>1
            dr=dr+getA(a-f)*[1;-0.04*(-1)^j];
        end
        
        line(pp(1,:) + dr(1), pp(2,:) + dr(2));    
    end
    
    f=getframe;

    % M(:,i) = getframe;    
    % imwrite(M(:,i),'img.png','PNG','Resolution',100);
    % mo= getframe;
    % mov=addframe(mov,mo);               
end

%%
% Energy

qabs=q*0;
v=zeros(size(q,1),2);
h=zeros(size(q,1),2);

T=0;
P=0;

for i=1:n
    v=zeros(size(q,1),2);
    h=zeros(size(q,1),2);
    for j=1:i
       qabs(:,i)=qabs(:,i)+q(:,j); 
       qabs(:,n+i)=qabs(:,n+i)+q(:,n+j);       
    end
    for j=1:i
       if j==i 
           b=-c{j,j};
       else
           b=-c{j,j}+c{j,j+1};
       end
       h=h+reshape(getA(qabs(:,j))*b,size(q,1),2);
       v=v-reshape(getA(qabs(:,j))*[0 1;-1 0]*b,size(q,1),2).*[qabs(:,n+j) qabs(:,n+j)];      
    end    
    T=T+mass{i}(1,1)*(v(:,1).^2+v(:,2).^2)*0.5+mass{i}(3,3)*qabs(:,n+i).^2*0.5;
    P=P+Model.mass{i}(1,1)*9.81*h(:,2);
end

figure; 
E0=T+P;
plot(t,T+P-E0(1));
title('Energy error');





##### SOURCE END #####
-->
   </body>
</html>